Introducción

El termino materia obscura hace referencia a todas las especies microbianas que no han sido caracterizadas. Estas especies impiden la caracterización de ecosistemas dominados por estos microorganismos, debido a la falta de entendimiento de sus funciones.

Redes tomando en cuenta los datos de todas las muestras.

Librerías utilizadas:

library(igraph)
library(igraphdata)
library(networkD3)
library(RCy3)
library(tidyverse)  
library(microbiomeDataSets)
library(mia)
library(phyloseq)
library(SpiecEasi)
library(ggplot2)
library(tidyr)

Datos a utilizar de “library(microbiomeDataSets)”:

datos_sprockett <- SprockettTHData()

ncol(datos_sprockett)
## [1] 2086

Proceso para construir el objeto phyloseq:

matriz_otu <- assays(datos_sprockett)$counts %>% as.matrix()
matriz_tax <- rowData(datos_sprockett) %>% as.data.frame() %>% as.matrix()
datos_muestras <- colData(datos_sprockett) %>% as.data.frame()

OTU <- otu_table(matriz_otu, taxa_are_rows = TRUE)
TAX <- tax_table(matriz_tax)
MUESTRAS <- sample_data(datos_muestras)
physeq <- phyloseq(OTU, TAX, MUESTRAS)

Separacion de objetos con y sin materia obscura:

physeq_especie <- tax_glom(physeq, taxrank = "Species", NArm = FALSE)

physeq_conocido <- subset_taxa(physeq_especie, !is.na(Species))  # sin materia obscura
physeq_todos <- physeq_especie  # con materia obscura

Filtrado por prevalencia (en este caso tomamos un filtrado del 10%, para obtener una red más densa):

filtrar_prevalencia <- function(physeq_objeto, umbral = 0.1) {
  prevalencia <- apply(otu_table(physeq_objeto), 1, function(x) mean(x > 0))
  conservar <- names(prevalencia[prevalencia >= umbral])
  prune_taxa(conservar, physeq_objeto)
}

physeq_conocido_filtrado <- filtrar_prevalencia(physeq_conocido, 0.1)
physeq_todos_filtrado    <- filtrar_prevalencia(physeq_todos, 0.1)

Renombramiento de especies para una red con mejor visualización e idenificación de NA’s:

renombrar_especies <- function(physeq_objeto) {
  tax <- tax_table(physeq_objeto)
  especies <- as.character(tax[, "Species"])
  
  especies_limpias <- ifelse(is.na(especies), "Desconocido", especies)
  
  nombres_unicos <- make.unique(especies_limpias)
  
  taxa_names(physeq_objeto) <- nombres_unicos
  return(physeq_objeto)
}

physeq_conocido_filtrado <- renombrar_especies(physeq_conocido_filtrado)
physeq_todos_filtrado    <- renombrar_especies(physeq_todos_filtrado)

Construcción de redes con SPIEC-EASI que utiliza el método de co-ocurrencia:

red_conocido <- spiec.easi(
  physeq_conocido_filtrado,
  method = "mb",
  lambda.min.ratio = 1e-1,
  nlambda = 20,
  sel.criterion = "bstars",
  pulsar.params = list(thresh = 0.1)
)

red_todos <- spiec.easi(
  physeq_todos_filtrado,
  method = "mb",
  lambda.min.ratio = 1e-1,
  nlambda = 20,
  sel.criterion = "bstars",
  pulsar.params = list(thresh = 0.1)
)

Transfromación a formato Igraph:

g_conocido <- adj2igraph(getRefit(red_conocido), vertex.attr = list(name = taxa_names(physeq_conocido_filtrado)))

g_todos    <- adj2igraph(getRefit(red_todos),    vertex.attr = list(name = taxa_names(physeq_todos_filtrado)))

Cálculo de métricas de centralidad:

calcular_metricas <- function(grafo) {
  list(
    nodos = vcount(grafo),
    aristas = ecount(grafo),
    grado_medio = mean(degree(grafo)),
    densidad = edge_density(grafo),
    clustering = transitivity(grafo, type = "global"),
    modularidad = modularity(cluster_louvain(grafo)),
    betweenness_media = mean(betweenness(grafo))
  )
}

metricas_conocido <- calcular_metricas(g_conocido)
metricas_todos    <- calcular_metricas(g_todos)

Tabla comparativa de metricas de centralidad:

tabla_comparativa <- tibble(
  Métrica              = names(metricas_todos),
  Con_materia_obscura  = unlist(metricas_todos),
  Sin_materia_obscura  = unlist(metricas_conocido)
)
print(tabla_comparativa)
## # A tibble: 7 × 3
##   Métrica           Con_materia_obscura Sin_materia_obscura
##   <chr>                           <dbl>               <dbl>
## 1 nodos                        115                  57     
## 2 aristas                      332                 124     
## 3 grado_medio                    5.77                4.35  
## 4 densidad                       0.0506              0.0777
## 5 clustering                     0.239               0.270 
## 6 modularidad                    0.516               0.533 
## 7 betweenness_media            101.                 45.6

Observación de redes:

plot(g_todos,
     vertex.size = degree(g_todos)*2,
     vertex.color = cluster_louvain(g_todos)$membership,
     vertex.label.cex = 0.7,
     edge.width = 1,
     main = "Red con Materia Obscura")

plot(g_conocido,
     vertex.size = degree(g_conocido)*2,
     vertex.color = cluster_louvain(g_conocido)$membership,
     vertex.label.cex = 0.7,
     edge.width = 1,
     main = "Red sin Materia Obscura")

Edición de las redes en Cytoscape para poder visualizarlas de mejor manera:

Red con materia obscura
Red con materia obscura
Red sin materia obscura
Red sin materia obscura

Red sin taxones obtenidos aleatoriamente

Creación de función para eliminar taxones de forma aleatoria:

eliminar_taxones_azar <- function(physeq, porcentaje = 0.3) {
  taxones <- taxa_names(physeq)
  eliminar <- sample(taxones, size = round(length(taxones) * porcentaje))
  prune_taxa(setdiff(taxones, eliminar), physeq)
}

physeq_reducido <- eliminar_taxones_azar(physeq_conocido_filtrado)

filtrar_prevalencia <- function(physeq_objeto, umbral = 0.1) {
  prevalencia <- apply(otu_table(physeq_objeto), 1, function(x) mean(x > 0))
  conservar <- names(prevalencia[prevalencia >= umbral])
  prune_taxa(conservar, physeq_objeto)
}
physeq_reducido_filt <- filtrar_prevalencia(physeq_reducido, 0.1)

Elaboración de red sin taxones al azar

se_reducido <- spiec.easi(
  physeq_reducido_filt,
  method = "mb",
  lambda.min.ratio = 1e-1,
  nlambda = 20,
  sel.criterion = "bstars",
  pulsar.params = list(thresh = 0.1)
)

g_reducido <- adj2igraph(getRefit(se_reducido), vertex.attr = list(name = taxa_names(physeq_reducido_filt)))


# Visualización de la red reducida
plot(g_reducido,
     vertex.size = degree(g_reducido)*2,
     vertex.color = cluster_louvain(g_reducido)$membership,
     vertex.label.cex = 0.7,
     edge.width = 1,
     main = "Red con 30% de taxones eliminados")

Calculo de métricas de la red

metricas <- function(g) {
  list(
    nodos = vcount(g),
    aristas = ecount(g),
    grado_medio = mean(degree(g)),
    densidad = edge_density(g),
    clustering = transitivity(g, type = "global"),
    modularidad = modularity(cluster_louvain(g))
  )
}

m_reducido <- metricas(g_reducido)

Edición de la red en Cytoscape para poder visualizarla de mejor manera:

Comparación de medidas de centralidad de todas las redes creadas hasta el momento

metricas_todos <- calcular_metricas(g_todos)
metricas_conocido <- calcular_metricas(g_conocido)
metricas_bootstrap <- calcular_metricas(g_reducido)  

tabla_comparativa_especie <- tibble(
  Métrica = names(metricas_todos),
  Con_materia_obscura = unlist(metricas_todos),
  Sin_materia_obscura = unlist(metricas_conocido),
  Bootstrap = unlist(metricas_bootstrap)
)

print(tabla_comparativa_especie)
## # A tibble: 7 × 4
##   Métrica           Con_materia_obscura Sin_materia_obscura Bootstrap
##   <chr>                           <dbl>               <dbl>     <dbl>
## 1 nodos                        115                  57        40     
## 2 aristas                      332                 124        54     
## 3 grado_medio                    5.77                4.35      2.7   
## 4 densidad                       0.0506              0.0777    0.0692
## 5 clustering                     0.239               0.270     0.247 
## 6 modularidad                    0.510               0.532     0.625 
## 7 betweenness_media            101.                 45.6      34.6

Visualización de gráfica comparativa

df_metricas_especie <- tabla_comparativa_especie %>%
  pivot_longer(cols = -Métrica, names_to = "Red", values_to = "Valor")

ggplot(df_metricas_especie, aes(x = Métrica, y = Valor, fill = Red)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +  #
  scale_fill_manual(values = c("Con_materia_obscura" = "lightblue", 
                               "Sin_materia_obscura" = "pink", 
                               "Bootstrap" = "purple")) +  # colores para cada red
  theme_minimal() +
  coord_flip()+
  labs(title = "Comparación de métricas entre redes",
       x = "Métrica",
       y = "Valor de métrica") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Redes con y sin materia obscura tomando en cuenta únicamente las muestras de heces

Filtrado de los datos.

physeq_feces <- subset_samples(physeq, Sample_Type == "Feces")

physeq_saliva <- subset_samples(physeq, Sample_Type == "Saliva") 

physeq_adultos <- subset_samples(physeq, Age_Class == "Adult")

physeq_infantes <- subset_samples(physeq, Age_Class == "Infant")

physeq_niños <- subset_samples(physeq, Age_Class == "Child")

Objeto Heces: Separar cada nuevo objeto con y sin material oscura.

physeq_feces <- subset_samples(physeq, Sample_Type == "Feces")


physeq_spp_f <- tax_glom(physeq_feces, taxrank = "Species", NArm = FALSE) 

physeq_sindm     <- subset_taxa(physeq_spp_f, !is.na(Species))  
physeq_conmd      <- physeq_spp_f

Filtrado por prevalencia:

#
prevalence_filter <- function(physeq_objeto, umbral = 0.1) {
  prevalencia <- apply(otu_table(physeq_objeto), 1, function(x) mean(x > 0))
  conservar <- names(prevalencia[prevalencia >= umbral])
  prune_taxa(conservar, physeq_objeto)
}
#
physeq_known_filt.feces <- prevalence_filter(physeq_sindm, 0.2) 
physeq_all_filt.feces   <- prevalence_filter(physeq_conmd, 0.2)

Renombramiento de especies para una red con mejor visualización e idenificación de NA’s:

#
renombrar_especies <- function(physeq_objeto) {
  tax <- tax_table(physeq_objeto)
  especies <- as.character(tax[, "Species"])
  
  # Reemplazar NA con "Desconocido"
  especies_limpias <- ifelse(is.na(especies), "Desconocido", especies)
  
  # Asegurar unicidad en los nombres
  nombres_unicos <- make.unique(especies_limpias)
  
  # Asignar los nuevos nombres como nombres de los taxones
  taxa_names(physeq_objeto) <- nombres_unicos
  return(physeq_objeto)
}
#
physeq_known_filt.feces <- renombrar_especies(physeq_known_filt.feces)
physeq_all_filt.feces    <- renombrar_especies(physeq_all_filt.feces)

Construcción de redes con SPIEC-EASI que utiliza el método de co-ocurrencia:

feces_known <- spiec.easi( 
  physeq_known_filt.feces,
  method = "mb",
  lambda.min.ratio = 1e-1,
  nlambda = 20,
  sel.criterion = "bstars",
  pulsar.params = list(thresh = 0.1)
)

feces_all <- spiec.easi(
  physeq_all_filt.feces,
  method = "mb",
  lambda.min.ratio = 1e-1,
  nlambda = 20,
  sel.criterion = "bstars",
  pulsar.params = list(thresh = 0.1)
)

Transformacion a formato igraph

feces_known <- adj2igraph(getRefit(feces_known), 
                          vertex.attr = list(name = taxa_names(physeq_known_filt.feces)))
feces_all   <- adj2igraph(getRefit(feces_all),   
                          vertex.attr = list(name = taxa_names(physeq_all_filt.feces)))

Cálculo de métricas de centralidad:

#
metricas <- function(g) {
  list(
    nodos = vcount(g),
    aristas = ecount(g),
    grado_medio = mean(degree(g)),
    densidad = edge_density(g),
    clustering = transitivity(g, type = "global"),
    modularidad = modularity(cluster_louvain(g))
  )
}

#
feces.m_known <- metricas(feces_known)
feces.m_all   <- metricas(feces_all)

Tabla comparativa de metricas de centralidad:

#COMPARACION:
tibble(
  Métrica        = names(feces.m_all),
  Con_materia_obscura = unlist(feces.m_all),
  Sin_materia_obscura = unlist(feces.m_known)
)
## # A tibble: 6 × 3
##   Métrica     Con_materia_obscura Sin_materia_obscura
##   <chr>                     <dbl>               <dbl>
## 1 nodos                   51                  27     
## 2 aristas                 92                  20     
## 3 grado_medio              3.61                1.48  
## 4 densidad                 0.0722              0.0570
## 5 clustering               0.237               0.140 
## 6 modularidad              0.551               0.579

Observación de redes:

# VISUALIZACIÓN con materia obscura)
plot(feces_all,
     vertex.size = degree(feces_all)*2,
     vertex.color = cluster_louvain(feces_all)$membership,
     vertex.label.cex = 0.7,
     edge.width = 1,
     main = "Red con Materia Obscura HECES")

# VISUALIZACIÓN DE g_known (sin materia obscura)
plot(feces_known,
     vertex.size = degree(feces_known)*2,
     vertex.color = cluster_louvain(feces_known)$membership,
     vertex.label.cex = 0.7,
     edge.width = 1,
     main = "Red sin Materia Obscura HECES")

Edición de las redes en Cytoscape para poder visualizarlas de mejor manera:

Red con materia oscura heces:

Red sin materia oscura heces:

Redes con y sin materia obscura tomando en cuenta únicamente las muestras de saliva

Objeto Saliva: Separar cada nuevo objeto con y sin material oscura.

physeq_saliva <- subset_samples(physeq, Sample_Type == "Saliva") 

physeq_spp_s <- tax_glom(physeq_saliva, taxrank = "Species", NArm = FALSE)

physeq_saliva.sin     <- subset_taxa(physeq_spp_s, !is.na(Species))  # sin materia oscura
physeq_saliva.com       <- physeq_spp_s  # con materia oscura

Filtrado por prevalencia:

#
prevalence_filter <- function(physeq_objeto, umbral = 0.1) {
  prevalencia <- apply(otu_table(physeq_objeto), 1, function(x) mean(x > 0))
  conservar <- names(prevalencia[prevalencia >= umbral])
  prune_taxa(conservar, physeq_objeto)
}
#
physeq_known_filt.saliva <- filtrar_prevalencia(physeq_saliva.sin, 0.2) #aplciar la funcion y sacar nuevos objetos
physeq_all_filt.saliva   <- filtrar_prevalencia(physeq_saliva.com, 0.2)

Renombramiento de especies para una red con mejor visualización e idenificación de NA’s:

#
renombrar_especies <- function(physeq_objeto) {
  tax <- tax_table(physeq_objeto)
  especies <- as.character(tax[, "Species"])
  
  # Reemplazar NA con "Desconocido"
  especies_limpias <- ifelse(is.na(especies), "Desconocido", especies)
  
  # Asegurar unicidad en los nombres
  nombres_unicos <- make.unique(especies_limpias)
  
  # Asignar los nuevos nombres como nombres de los taxones
  taxa_names(physeq_objeto) <- nombres_unicos
  return(physeq_objeto)
}
#
physeq_known_filt.saliva <- renombrar_especies(physeq_known_filt.saliva)
physeq_all_filt.saliva    <- renombrar_especies(physeq_all_filt.saliva)

Construcción de redes con SPIEC-EASI que utiliza el método de co-ocurrencia:

# REDES SPIEC-EASI
#analsiis de redes de co-ocurrencia 
saliva_known <- spiec.easi( #NO DM
  physeq_known_filt.saliva,
  method = "mb",
  lambda.min.ratio = 1e-1,
  nlambda = 20,
  sel.criterion = "bstars",
  pulsar.params = list(thresh = 0.1)
)
## Applying data transformations...
## Selecting model with pulsar using bstars...
## Fitting final estimate with mb...
## done
saliva_all <- spiec.easi(
  physeq_all_filt.saliva,
  method = "mb",
  lambda.min.ratio = 1e-1,
  nlambda = 20,
  sel.criterion = "bstars",
  pulsar.params = list(thresh = 0.1)
)
## Applying data transformations...
## Selecting model with pulsar using bstars...
## Fitting final estimate with mb...
## done

Transformacion a formato igraph

saliva_known <- adj2igraph(getRefit(saliva_known), 
                           vertex.attr = list(name = taxa_names(physeq_known_filt.saliva)))
saliva_all   <- adj2igraph(getRefit(saliva_all),   
                           vertex.attr = list(name = taxa_names(physeq_all_filt.saliva)))

Cálculo de métricas de centralidad:

metricas <- function(g) {
  list(
    nodos = vcount(g),
    aristas = ecount(g),
    grado_medio = mean(degree(g)),
    densidad = edge_density(g),
    clustering = transitivity(g, type = "global"),
    modularidad = modularity(cluster_louvain(g))
  )
}
#
saliva.m_known <- metricas(saliva_known)
saliva.m_all   <- metricas(saliva_all)

Tabla comparativa de metricas de centralidad:

#COMPARACION:
tibble(
  Métrica        = names(saliva.m_all),
  Con_materia_obscura = unlist(saliva.m_all),
  Sin_materia_obscura = unlist(saliva.m_known)
)
## # A tibble: 6 × 3
##   Métrica     Con_materia_obscura Sin_materia_obscura
##   <chr>                     <dbl>               <dbl>
## 1 nodos                   75                  54     
## 2 aristas                204                 112     
## 3 grado_medio              5.44                4.15  
## 4 densidad                 0.0735              0.0783
## 5 clustering               0.221               0.265 
## 6 modularidad              0.484               0.520

Observación de redes:

plot(saliva_all,
     vertex.size = degree(saliva_all)*2,
     vertex.color = cluster_louvain(saliva_all)$membership,
     vertex.label.cex = 0.7,
     edge.width = 1,
     main = "Red con Materia Obscura saliva")

# VISUALIZACIÓN DE g_known (sin materia obscura)
plot(saliva_known,
     vertex.size = degree(saliva_known)*2,
     vertex.color = cluster_louvain(saliva_known)$membership,
     vertex.label.cex = 0.7,
     edge.width = 1,
     main = "Red sin Materia Obscura saliva")

Edición de las redes en Cytoscape para poder visualizarlas de mejor manera:

Red con materia oscura saliva:

Red sin materia oscura saliva:

## Redes con y sin materia obscura tomando en cuenta únicamente las muestras de niños

Objeto niños: Separar cada nuevo objeto con y sin material oscura.

physeq_niños <- subset_samples(physeq, Age_Class == "Child")

physeq_spp_niños <- tax_glom(physeq_niños, taxrank = "Species", NArm = FALSE)

physeq_niños.sin     <- subset_taxa(physeq_spp_niños, !is.na(Species))  # sin materia oscura
physeq_niños.com       <- physeq_spp_niños  # con materia oscura

Filtrado por prevalencia:

#
prevalence_filter <- function(physeq_objeto, umbral = 0.1) {
  prevalencia <- apply(otu_table(physeq_objeto), 1, function(x) mean(x > 0))
  conservar <- names(prevalencia[prevalencia >= umbral])
  prune_taxa(conservar, physeq_objeto)
}
#
physeq_conocido_filtrado_niños <- filtrar_prevalencia(physeq_niños.sin, 0.1)
physeq_todos_filtrado_niños    <- filtrar_prevalencia(physeq_niños.com, 0.1)

Renombramiento de especies para una red con mejor visualización e idenificación de NA’s:

#
renombrar_especies <- function(physeq_objeto) {
  tax <- tax_table(physeq_objeto)
  especies <- as.character(tax[, "Species"])
  
  # Reemplazar NA con "Desconocido"
  especies_limpias <- ifelse(is.na(especies), "Desconocido", especies)
  
  # Asegurar unicidad en los nombres
  nombres_unicos <- make.unique(especies_limpias)
  
  # Asignar los nuevos nombres como nombres de los taxones
  taxa_names(physeq_objeto) <- nombres_unicos
  return(physeq_objeto)
}
#

# Renombrar los taxa en physeq_especie (aplicable a physeq_todos y physeq_conocido)
physeq_conocido_filtrado_niños <- renombrar_especies(physeq_conocido_filtrado_niños)
physeq_todos_filtrado_niños    <- renombrar_especies(physeq_todos_filtrado_niños)

Construcción de redes con SPIEC-EASI que utiliza el método de co-ocurrencia:

# CONSTRUCCIÓN DE REDES CON SPIEC-EASI
red_conocido_niños <- spiec.easi(
  physeq_conocido_filtrado_niños,
  method = "mb",
  lambda.min.ratio = 1e-1,
  nlambda = 20,
  sel.criterion = "bstars",
  pulsar.params = list(thresh = 0.1)
)

red_todos_niños <- spiec.easi(
  physeq_todos_filtrado_niños,
  method = "mb",
  lambda.min.ratio = 1e-1,
  nlambda = 20,
  sel.criterion = "bstars",
  pulsar.params = list(thresh = 0.1)
)

Transformacion a formato igraph.

g_conocido.niños <- adj2igraph(getRefit(red_conocido_niños), vertex.attr = list(name = taxa_names(physeq_conocido_filtrado_niños)))
g_todos.niños    <- adj2igraph(getRefit(red_todos_niños),    vertex.attr = list(name = taxa_names(physeq_todos_filtrado_niños)))

Cálculo de métricas de centralidad:

metricas <- function(g) {
  list(
    nodos = vcount(g),
    aristas = ecount(g),
    grado_medio = mean(degree(g)),
    densidad = edge_density(g),
    clustering = transitivity(g, type = "global"),
    modularidad = modularity(cluster_louvain(g))
  )
}
#
metricas_conocido_niños <- calcular_metricas(g_conocido.niños)
metricas_todos_niños    <- calcular_metricas(g_todos.niños)

Tabla comparativa de metricas de centralidad:

tabla_comparativa.niños <- tibble(
  Métrica              = names(metricas_todos_niños),
  Con_materia_obscura  = unlist(metricas_todos_niños),
  Sin_materia_obscura  = unlist(metricas_conocido_niños)
)
print(tabla_comparativa.niños)
## # A tibble: 7 × 3
##   Métrica           Con_materia_obscura Sin_materia_obscura
##   <chr>                           <dbl>               <dbl>
## 1 nodos                        162                  62     
## 2 aristas                      531                  73     
## 3 grado_medio                    6.56                2.35  
## 4 densidad                       0.0407              0.0386
## 5 clustering                     0.102               0.149 
## 6 modularidad                    0.414               0.617 
## 7 betweenness_media            161.                 55.1

Observación de redes:

plot(g_todos.niños,
     vertex.size = degree(g_todos.niños)*2,
     vertex.color = cluster_louvain(g_todos.niños)$membership,
     vertex.label.cex = 0.7,
     edge.width = 1,
     main = "Red con Materia Obscura en niños")

plot(g_conocido.niños,
     vertex.size = degree(g_conocido.niños)*2,
     vertex.color = cluster_louvain(g_conocido.niños)$membership,
     vertex.label.cex = 0.7,
     edge.width = 1,
     main = "Red sin Materia Obscura en niños")

Edición de las redes en Cytoscape para poder visualizarlas de mejor manera:

Con materia oscura en niños (heces y saliva):

Sin materia oscura en niños (heces y saliva):

Redes con y sin materia obscura tomando en cuenta únicamente las muestras de infantes

Objeto infantes: Separar cada nuevo objeto con y sin material oscura.

physeq_infantes <- subset_samples(physeq, Age_Class == "Infant")

physeq_spp_infantes <- tax_glom(physeq_infantes, taxrank = "Species", NArm = FALSE)

physeq_infantes.sin     <- subset_taxa(physeq_spp_infantes, !is.na(Species))  # sin materia oscura
physeq_infante.com       <- physeq_spp_infantes  # con materia oscura

Filtrado por prevalencia:

#
prevalence_filter <- function(physeq_objeto, umbral = 0.1) {
  prevalencia <- apply(otu_table(physeq_objeto), 1, function(x) mean(x > 0))
  conservar <- names(prevalencia[prevalencia >= umbral])
  prune_taxa(conservar, physeq_objeto)
}
#
physeq_conocido_filtrado.infantes <- filtrar_prevalencia(physeq_infantes.sin, 0.1)
physeq_todos_filtrado.infantes   <- filtrar_prevalencia(physeq_infante.com, 0.1)

Renombramiento de especies para una red con mejor visualización e idenificación de NA’s:

#
renombrar_especies <- function(physeq_objeto) {
  tax <- tax_table(physeq_objeto)
  especies <- as.character(tax[, "Species"])
  
  # Reemplazar NA con "Desconocido"
  especies_limpias <- ifelse(is.na(especies), "Desconocido", especies)
  
  # Asegurar unicidad en los nombres
  nombres_unicos <- make.unique(especies_limpias)
  
  # Asignar los nuevos nombres como nombres de los taxones
  taxa_names(physeq_objeto) <- nombres_unicos
  return(physeq_objeto)
}
#
# Renombrar los taxa en physeq_especie (aplicable a physeq_todos y physeq_conocido)
physeq_conocido_filtrado.infantes <- renombrar_especies(physeq_conocido_filtrado.infantes)
physeq_todos_filtrado.infantes    <- renombrar_especies(physeq_todos_filtrado.infantes)

Construcción de redes con SPIEC-EASI que utiliza el método de co-ocurrencia:

# CONSTRUCCIÓN DE REDES CON SPIEC-EASI
red_conocido_infantes <- spiec.easi(
  physeq_conocido_filtrado.infantes,
  method = "mb",
  lambda.min.ratio = 1e-1,
  nlambda = 20,
  sel.criterion = "bstars",
  pulsar.params = list(thresh = 0.1)
)

red_todos_infantes <- spiec.easi(
  physeq_todos_filtrado.infantes,
  method = "mb",
  lambda.min.ratio = 1e-1,
  nlambda = 20,
  sel.criterion = "bstars",
  pulsar.params = list(thresh = 0.1)
)

Transformacion a formato igraph.

g_conocido.infantes <- adj2igraph(getRefit(red_conocido_infantes), vertex.attr = list(name = taxa_names(physeq_conocido_filtrado.infantes)))

g_todos.infantes  <- adj2igraph(getRefit(red_todos_infantes),    vertex.attr = list(name = taxa_names(physeq_todos_filtrado.infantes)))

Cálculo de métricas de centralidad:

metricas <- function(g) {
  list(
    nodos = vcount(g),
    aristas = ecount(g),
    grado_medio = mean(degree(g)),
    densidad = edge_density(g),
    clustering = transitivity(g, type = "global"),
    modularidad = modularity(cluster_louvain(g))
  )
}
#
metricas_conocido_infantes <- calcular_metricas(g_conocido.infantes)
metricas_todos_infantes    <- calcular_metricas(g_todos.infantes)

Tabla comparativa de metricas de centralidad:

tabla_comparativa.infantes <- tibble(
  Métrica              = names(metricas_todos_infantes),
  Con_materia_obscura  = unlist(metricas_todos_infantes),
  Sin_materia_obscura  = unlist(metricas_conocido_infantes)
)
print(tabla_comparativa.infantes)
## # A tibble: 7 × 3
##   Métrica           Con_materia_obscura Sin_materia_obscura
##   <chr>                           <dbl>               <dbl>
## 1 nodos                         87                  54     
## 2 aristas                      218                  74     
## 3 grado_medio                    5.01                2.74  
## 4 densidad                       0.0583              0.0517
## 5 clustering                     0.245               0.214 
## 6 modularidad                    0.504               0.518 
## 7 betweenness_media             62.0                40.0

Observación de redes:

plot(g_todos.infantes,
     vertex.size = degree(g_todos.infantes)*2,
     vertex.color = cluster_louvain(g_todos.infantes)$membership,
     vertex.label.cex = 0.7,
     edge.width = 1,
     main = "Red con Materia Obscura en infantes")

plot(g_conocido.infantes,
     vertex.size = degree(g_conocido.infantes)*2,
     vertex.color = cluster_louvain(g_conocido.infantes)$membership,
     vertex.label.cex = 0.7,
     edge.width = 1,
     main = "Red sin Materia Obscura en infantes")

Edición de las redes en Cytoscape para poder visualizarlas de mejor manera:

Con materia oscura en infantes (heces y saliva):

Sin materia oscura infantes (saliva y heces):

Redes con y sin materia obscura tomando en cuenta únicamente las muestras de adultos

Objeto Adultos: Separar cada nuevo objeto con y sin material oscura.

physeq_adultos <- subset_samples(physeq, Age_Class == "Adult")

physeq_spp_adultos <- tax_glom(physeq_adultos, taxrank = "Species", NArm = FALSE)

physeq_adultos.sin     <- subset_taxa(physeq_spp_adultos, !is.na(Species))  # sin materia oscura
physeq_adultos.com       <- physeq_spp_adultos  # con materia oscura

Filtrado por prevalencia:

#
prevalence_filter <- function(physeq_objeto, umbral = 0.1) {
  prevalencia <- apply(otu_table(physeq_objeto), 1, function(x) mean(x > 0))
  conservar <- names(prevalencia[prevalencia >= umbral])
  prune_taxa(conservar, physeq_objeto)
}
#
physeq_conocido_filtrado.adultos <- filtrar_prevalencia(physeq_adultos.sin, 0.1)
physeq_todos_filtrado.adultos   <- filtrar_prevalencia(physeq_adultos.com, 0.1)

Renombramiento de especies para una red con mejor visualización e idenificación de NA’s:

#
renombrar_especies <- function(physeq_objeto) {
  tax <- tax_table(physeq_objeto)
  especies <- as.character(tax[, "Species"])
  
  # Reemplazar NA con "Desconocido"
  especies_limpias <- ifelse(is.na(especies), "Desconocido", especies)
  
  # Asegurar unicidad en los nombres
  nombres_unicos <- make.unique(especies_limpias)
  
  # Asignar los nuevos nombres como nombres de los taxones
  taxa_names(physeq_objeto) <- nombres_unicos
  return(physeq_objeto)
}
#
# Renombrar los taxa en physeq_especie (aplicable a physeq_todos y physeq_conocido)
physeq_conocido_filtrado.adultos <- renombrar_especies(physeq_conocido_filtrado.adultos)
physeq_todos_filtrado.adultos    <- renombrar_especies(physeq_todos_filtrado.adultos)

Construcción de redes con SPIEC-EASI que utiliza el método de co-ocurrencia:

# CONSTRUCCIÓN DE REDES CON SPIEC-EASI
red_conocido_adultos <- spiec.easi(
  physeq_conocido_filtrado.adultos,
  method = "mb",
  lambda.min.ratio = 1e-1,
  nlambda = 20,
  sel.criterion = "bstars",
  pulsar.params = list(thresh = 0.1)
)
## Applying data transformations...
## Selecting model with pulsar using bstars...
## Fitting final estimate with mb...
## done
red_todos_adultos <- spiec.easi(
  physeq_todos_filtrado.adultos,
  method = "mb",
  lambda.min.ratio = 1e-1,
  nlambda = 20,
  sel.criterion = "bstars",
  pulsar.params = list(thresh = 0.1)
)
## Applying data transformations...
## Selecting model with pulsar using bstars...
## Fitting final estimate with mb...
## done

Transformacion a formato igraph.

g_conocido.adultos <- adj2igraph(getRefit(red_conocido_adultos), 
                                 vertex.attr = list(name = taxa_names(physeq_conocido_filtrado.adultos)))

g_todos.adultos  <- adj2igraph(getRefit(red_todos_adultos),    
                               vertex.attr = list(name = taxa_names(physeq_todos_filtrado.adultos)))

Cálculo de métricas de centralidad:

metricas <- function(g) {
  list(
    nodos = vcount(g),
    aristas = ecount(g),
    grado_medio = mean(degree(g)),
    densidad = edge_density(g),
    clustering = transitivity(g, type = "global"),
    modularidad = modularity(cluster_louvain(g))
  )
}
#
metricas_conocido_adultos <- calcular_metricas(g_conocido.adultos)
metricas_todos_adultos    <- calcular_metricas(g_todos.adultos)

Tabla comparativa de metricas de centralidad:

tabla_comparativa.adultos <- tibble(
  Métrica              = names(metricas_todos_adultos),
  Con_materia_obscura  = unlist(metricas_todos_adultos),
  Sin_materia_obscura  = unlist(metricas_conocido_adultos)
)
print(tabla_comparativa.adultos)
## # A tibble: 7 × 3
##   Métrica           Con_materia_obscura Sin_materia_obscura
##   <chr>                           <dbl>               <dbl>
## 1 nodos                        177                  73     
## 2 aristas                      737                 174     
## 3 grado_medio                    8.33                4.77  
## 4 densidad                       0.0473              0.0662
## 5 clustering                     0.202               0.262 
## 6 modularidad                    0.446               0.428 
## 7 betweenness_media            164.                 66.6

Observación de redes:

plot(g_todos.adultos,
     vertex.size = degree(g_todos.adultos)*2,
     vertex.color = cluster_louvain(g_todos.adultos)$membership,
     vertex.label.cex = 0.7,
     edge.width = 1,
     main = "Red con Materia Obscura en adultos")

plot(g_conocido.adultos,
     vertex.size = degree(g_conocido.adultos)*2,
     vertex.color = cluster_louvain(g_conocido.adultos)$membership,
     vertex.label.cex = 0.7,
     edge.width = 1,
     main = "Red sin Materia Obscura en infantes")

Edición de las redes en Cytoscape para poder visualizarlas de mejor manera:

Con materia oscura en adultos (heces y saliva):

Sin materia oscura adultos (saliva y heces):